function Vmat = European_put_imp(r,sigma,T,Smax,K)
%  Discussion:
%
%    This program solves black-scholes equation using FD 
%    with an implicit scheme.
%
%    Inputs: r--interest rate, sigma--volatility, T--Maturity,Smax--max
%    underlying price, K--strike
%    Output: Option Price of European Put Option
%%%%%%%%%Input Parameters%%%%%%%%%
if (nargin<1)
    r = 0.02;
    sigma = 0.4;
    T = 1;
    Smax = 10;
    K = 5;
else
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 101;
M = 31;
dt = T/(N-1);
ds = Smax/(M-1);
k = dt/ds^2;
l = dt/ds;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
v = zeros(M,1);
Vmat = zeros ( M, N );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S = linspace(0,Smax,M);
a = -0.5*(sigma^2*S.^2*k+r*S*l);
b = 1 + sigma^2*S.^2*k + r*dt;
c = -0.5*(sigma^2*S.^2*k-r*S*l);
  for j = 1 : N

    if ( j == 1 )
      v(:) = max(K-S',0);
    else
      A = diag([0,a(1:M-2)],1)+diag([1,b(1:M-2),1])+diag([c(1:M-2),0],-1);
      v = inv(A)*v;
    end
    v(M)= 0;
    v(1) = K*exp(-r*(j-1)*dt);
    Vmat(:,j) = v;
  end
  % Visualization
[X,Y] = meshgrid(0:dt:1,S);
surf(X,Y,Vmat)
end